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We use path-integral Monte Carlo (PIMC) to study the effects of adding a long-range 
repulsive Coulomb interaction to the usual Van der Waals interaction between two atoms 
of a submonolayer quantum film such as helium on graphite or a pure two-dimensional 
superfluid. Such interactions frustrate or compete with the natural tendency of the sys- 
tem for phase separation namely to form a macroscopic liquid or solid phase. We find 
that as a function of the relative strength of the long-range repulsion, surface coverage 
and temperature, the system undergoes a series of transformations, including a trian- 
gular Wigner-like crystal of clusters, a charge stripe-ordered phase and a fiuid phase. 
The goal of these studies is to understand the role of quantum fluctuations when such 
competing interactions appear together with formation of preexisting electron pairs as 
might be the case in cuprate superconductors. 
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1. INTRODUCTION 

In a number of physical systems where there are competing interactions such as a 
short-range attractive part and a long-range repulsion, a rich phase diagram is ex- 
pected. In the particular case where the long-range repulsion is the usual Coulomb 
interaction, nuclear matter and electrons in a solid such as the cuprate supercon- 
ductors are examples of this case. In the latter case, the electrons in the material 
can find themselves experiencing an attractive interaction which might also lead to 
pairing or phase separation, but when they form clusters the accumulated charge 
cannot exceed a threshold value due to the Coulomb repulsion. 

In this work we focus our attention to the case of a two dimensional quantum 
fluid of charged bosons with short-range attractive interaction and a hard core re- 
pulsion. This model, for example, could be used to understand certain aspects of 
the conduction electrons of the copper-oxide plane of the cuprate superconductors. 
It has been found^'^'^ that in cuprate superconductors, as modeled by the so-called 
t — J model, holes tend to phase separate in a hole rich and an antiferromagnetically 
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ordered region. A pair of holes in a quantum antiferromagnet is strongly bound^ 
and in fact the hole rich phase can be understood as a fluid of weakly attracted 
bosons- the bound pairs of holcs^. In this paper, wc consider such a system of bosons 
which interact via a short range attractive interaction, a hard core repulsion and 
a long-range Coulomb repulsion. A system of interacting bosons is free from the 
infamous minus problem which hinders the quantum simulation of fermions. There- 
fore we should be able to address the phase diagram of this model without serious 
difficulties. 

In the present paper our goal is to understand the influence of a weak Coulomb 
repulsion on a quantum Van dcr Waals-Bosc liquid monolayer film. Our model of 
the neutral Van der Waals-Bose liquid is a liquid helium film, because we have 
gained understanding of the phases of this system from experimental studies®'^'^ 
and from our previous PIMC studies^. More importantly, using the PIMC method 
most of the observed features of these films and the phase diagram of each layer as a 
function of surface coverage can be reproduced with detail. In these PIMC studies ^ 
of submonolayer, second, third and higher layers of helium on graphite well-known 
helium-helium potential and helium-graphite interaction were used without the in- 
troduction of any fitting parameter. In this paper, we begin from the Hamiltonian 
that describes a submonolayer helium film on graphite and we modify it by adding 
a weak Coulomb repulsion term to the helium-helium interaction potential which 
depends on a coupling constant a. Wc will use PIMC method to study such a film 
as a function of temperature, surface density and the parameter a. 

2. THE APPROACH 

For particle-substrate interaction we used the laterally averaged helium-graphite 
interaction^^ used in previous PIMC studies of helium adsorption on graphite^. In 
order to study the role of a long-range Coulomb repulsion on the Van der Waals 
helium- helium atom interaction we added a 1/r term to a Lennard- Jones interaction 
potential as follows: 

V{r) = -4e((^)' - i^)'') + aVc{L,r) (1) 

where Vc{L,r) is the 1/r modified to include the influence of the periodic image 
charges in a finite-size L x L lattice with periodic boundary conditions. This is 

implemented for two-dimensions by the well-known Ewald summation mcthod^^. 
For any given value of the parameter a we adjust the value of e in order for the 
minimum of the Lennard-Jones potential to be at the same value Vmin = W°K. 
This is the approximate minimum value of the Aziz potential"^''^ that describes the 
helium- helium interaction and it is Tised to study neutral helium quantum films. The 
total interaction potential is shown in Fig. 1 for various values of the parameter a. 
Notice that in order to keep Vmin a constant, increasing the value of a requires an 
increased value of e and this leads to a harder core potential and a more attractive 
force dragging the neighbor to toward the potential minimum. 



2, 2008 17:49 WSPC/INSTRUCTION FILE paper 



Monolayer charged quantum films 3 

We will use the path-integral Monte Carlo method (PIMC) for superfluid films 
as developed in Ref.^. We have carried out simulations at T = 0.5°.^' and T = 1°K 

for all 100 combinations of 10 density values starting from p = 0.0120A~^ up to 
0.04801-2 with 10 values of a starting from a = 50°KA up to a = 500°KA in 
increments of 50°KA. We kept the number of particles fixed at 30 and we varied 
the size of the simulation cell in order to vary the density. We have computed the 
pair distribution function, and contour plots of the particle probability density for 
all of the above parameter values. Approximately 10,000 Monte Carlo sweeps were 
used for thermalization and about 10,000 sweeps for measurements of observable. 
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Fig. 1. The total interaction potential is shown for various values of the parameter a. The value 
of e is adjusted such that the minimum of the interaction is fixed at —10°K. 



3. RESULTS 

In Fig. 2 we give a thermalized space-time configuration of the system, contour 

plots of the probability distribution averaged over 10,000 such configurations and 
the pair distribution function for surface density p = 0.01201"^ and temperature 
T = 500mK for various values of a. We notice that while for this surface density 
in the absence of the Coidomb repulsion the system phase separates, the presence 
of even weak Coulomb interaction tends to segregate the particles into clusters 
which spread over the entire film surface. In addition, the pair distribution function 
develops a second peak at the average inter-cluster distance. There is still a peak 
near the minimum of the Van der Waals attractive part of the interaction but for 
small a this peak is much smaller compared to that of the neutral fiuid. Notice that 
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Fig. 2. A thermalizcd space-time configuration, the contour plot and the pair distribution at 
surface density 0.0120A-2, T = 0.5°K and for a = (top row) a = 50° XA (second row from the 
top) a = 100° X A (third from top row) and a = 350° X A (fourth row from the top) 
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as we raise the strength of the Coulomb repulsion the system forms a triangular 
superlattice of clusters. The pair distribution function has a sharp peak near the 

minimum of the Van der Waals attraction and an oscillatory behavior with a period 
equal to the distance between nearest neighbor clusters. 
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Fig. 3. A thermalized space-time configuration, the contour plot and the pair distribution at 
surface density O.OI20A-2, T = 1°K and a = 50°KA (top row) and a = 35Q°KA (bottom row) 

By raising the temperature to 1°K the picture remains similar as shown in 
Fig. 3 with the tendency for the clusters and the particles to fluctuate more which 
will lead to the melting of the clusters at higher temperature. It is not clear if the 
system will become an ordered superlattice of such clusters or will remain a fluid 
at zero temperature at small values of a. It is well known that, strictly speaking, 
a two-dimensional system such as ours cannot break a continuous symmetry and 
crystallize at any finite temperature; however, a finite-size system may appear as 
ordered because the correlation length for low enough temperature becomes larger 
than the system size. We have observed that, as a function of the simulation time, 
a finite-size 2D system drifts as a whole. As a result, while within a relatively short 
time-scale the system appears ordered if we average configurations over a much 
larger time scale, the system will appear disordered in a contour plot. This time 
scale depends on the value of a. For small values of a it is easier for the system 
to migrate as a whole while keeping the particle correlations. This is demonstrated 
in the figures for a = 50, 100°KA where, while the particle correlations show that 
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Fig. 4. A thcrmalized space-time configuration, the contour plot and the pair distribution at 
surface density 0.0245A~^, T = bOOmK and a = (top row), a = bO°KA (second from the top 
row), a = 100° K A and a = 300° /C A (bottom row) 
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Fig. 5. A thcrmalizcd space-time configuration, the contour plot and the pair distribution at 
surface density 0.04801-^, T = 500mX and a = (top row), a = 200° K A (second from the top 
row), a = 300° if A and a = 350° K A (botto m row) 
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Fig. 6. A thermalized space-time configuration, the contour plot and the pair distribution at 
surface density 0.0480A-2, a = 450° K A and T = 500mK (top row), T = lOOOmK (bottom row) 



they are those of a Wigner crystal of clusters, the contour plots, taken over 10,000 
configurations, indicate that the system is uniform. The larger values of a require 
averaging over much larger number of configurations to see the effects of the drift 
discussed previously. 

In Fig. 4 a thermalized snapshot of a space-time configuration along with a 
contour plot of the probability distribution and the pair distribution is shown for 
a higher density of 0.0245A~^, T = 500mK and various values of the Coulomb 
repulsion parameter a. The situation is similar to that of density p = 0.0120A~^. 
We note, however, that while the inter-cluster distance did not change significantly 
the number of atoms in each cluster has increased by approximately the factor of 
the density increase. 

In Fig. 5 we present the same three functions for p = 0.0480A~^ which is very 
close to the equilibrium density of the neutral liquid. Low values of the strength of 
Coulomb repulsion term have small influence on these functions including the pair 
distribution function. We need to increase to a value a = 200° K A in order to see an 
influence of the repulsion on g(r) comparable to that occurring at a = 50° K A for 
lower values of density. The reason for this insensitivity is the fact at the equilibrium 
density the compressibility of the fluid is finite, while that of the phase separated 
state is close to zero. 

Notice as the values of a increases the clustering that occurs is clusters of vacan- 
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Fig. 7. The phase diagram of a submonolayer Van der Waals-Bose quantum film where an addi- 
tional long-range Coulomb interaction has been added. 



cies which tend to order on a triangular superlattice. We have not clearly observed 
an ordered superlattice of clusters of vacancies. This tendency for ordering occurs 
near the equilibrium density where in a phase separated state there is much more 
fluid than vacant space. The Coulomb repulsion tends to compress and spread the 
fluid equally in all directions thus creating clusters of particles or of vacancies or 
stripes. The stripes occiir at densities near the equilibrium density and for large 
values of a. Notice that g(j) has three peaks when the stripes form or tend to 
form in Fig. 5. The short-distance peak is that which corresponds to the Van der 
Waals minimum, the peak at the longest distance corresponds to the inter-stripe 
distance and the intermediate peak corresponds to the a further ordering within the 
stripe. Notice in the contour plots of Fig. 5 that the stripes are formed by merging 
of clusters and the stripe width is modulated along the stripe and the modulation 
periodicity defines this third peak structure of the pair distribution fmiction. 

In Fig. 6 we compare the situation for a larger value of a = 450°/'irA for the equi- 
librium density and at temperature T = SOOmif and T = lOOOmfc. Notice that the 
detailed ordering of particles with a stripe begins to go away at the higher temper- 
ature because of thermal fluctuations. However, the peaks of the pair distribution 
function remain unchanged. 
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4. PHASE DIAGRAM-CONCLUSIONS 

Our results are summarized on the phase diagram shown in Fig. 7. Within our large 
error bars in the determination of the phase boundaries of the three phases, we have 

not been able to find any significant shift in the phase boundaries with temperature. 
The phase boundaries separate a uniform fluid phase from a phase of bubbles or 
clusters ordered in a triangular lattice formation and a phase of ordered stripes at 
high surface density and relatively strong Coulomb repulsion strength. 

In order to draw conclusions relevant to the case of cuprate superconductors we 
need to also study the superfluid density of the system as a function of the density 
and the Coulomb repulsion. Due to the non-local nature of the permutation cycles 
and the winding number operator this requires a significantly larger computational 
time scales than the calculation presented here and we are in the process of carrying 
it out. 
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